/*  scomdiv.cpp  by K.Tsuru  */
// function ID = 904
/**************************
SComplex class
division (*this)/z = (a+bi)/(c+di)
 (a+bi)(c-di)   ac+bd +(bc-ad)i  a*(c/d)+b + {b*(c/d)-a}i    a*r+b + (b*r-a)i
=------------ = ------------- = ------------------------- = ----------------- (r = c/d)
  c^2+d^2       z.Norm()             c*(c/d) + d                 c*r + d

On May 16, 2010
Tested again
On June 20, 2015
Changed onece more by the
[reference] Kouya, Tomonori "Shohokarano FFT" chap. 3.
url : http://na-inet.jp/fft/ (in Japanese)
**************************/
#ifndef SN_H
#include "sn.h"
#endif

static const char* func = "SComplex /";

SComplex& SComplex::operator/=(const SComplex& z) 
{
    if ( z.IsZero(904) ) im.SetError(im.DIVIDED_BY_ZERO, func, 904);
    // (a+bi)/c
    if (z.im.Sign(904) == SNumber::ZERO) return ((*this) /= z.re); // "z" is real number.

    // (a+bi)/ci = b/c -(a/c)i
    if (z.re.Sign(904) == SNumber::ZERO) {
        const SDouble r = im / z.im;
        im = -re / z.im;
        re = r;
        return *this;
    }
#if 0
    const SDouble r = re * z.re + im * z.im;
    const SDouble recNorm_z = DReciprocal( z.Norm() );  // = 1/(c^2 + d^2)
    im = ( im * z.re - re * z.im ) * recNorm_z;
    re = r * recNorm_z;
#else
    SDouble r, s, w;
    if(z.re >= z.im) {
        r  = z.im / z.re; s = z.re + z.im * r;
        w  = (re + im * r) / s;
        im = (-re * r + im) / s;
        re = w;
    }else {
        r = z.re / z.im; s = z.re * r + z.im;
        w  = (re * r + im) / s;
        im = (-re + im * r)/ s;
        re = w;
    }
#endif
    return *this;
}